%matplotlib inline
import os,sys
sys.path.append('/home/dpirvu/python_stuff/')
sys.path.append('./bubbles_codes/')
from plotting import *
from bubble_tools import *
from experiment import *
from celluloid import Camera
bubbleList, velocitesList, instantonList, tmpList, fldcritList, tcritList, encritList = [], [], [], [], [], [], []
for tmp in [0,1,2,3]:
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
bubbleList.append(np.load(average_file(*exp_params)))
velocitesList.append(np.load(velocities_file(*exp_params)))
instantonList.append(np.load(supercrit_instanton_file(*exp_params)))
tmpList.append(tmp)
fldcritList.append(np.load(critfield_file(*exp_params)))
tcritList.append(np.load(crittimes_file(*exp_params)))
encritList.append(np.load(critenerg_file(*exp_params)))
# standardize plots
lsl = lambda tmp: ('-' if tmp%4==0 else '-.' if tmp%4==1 else '--' if tmp%4==3 else ':')
free_eigenbasis_mom = lambda la, ph0: np.asarray([norm(ph0)*(w2(la)[k]**0.25) if kk!=0. else 0. for k,kk in enumerate(klist)])
thermal_eigenbasis_mom = lambda la, ph0, te: free_eigenbasis_mom(la, ph0) * np.sqrt(2./(np.exp(w2(la)**0.5/te)-1.))
fluct_stdev_mom = lambda la, ph0, te: np.sqrt(np.sum(np.abs(thermal_eigenbasis_mom(la, ph0, te))**2.))
%run './bubbles_codes/plotting.py'
%run './bubbles_codes/bubble_tools.py'
if True:
tmp = 0
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
sigmamom = fluct_stdev_mom(lamb, phi0, temp)
print(sigmafld, sigmamom, right_Vmax.x, m2(lamb)**0.5, temp/m2(lamb)**0.5)
print((right_Vmax.x-np.pi)/sigmafld)
philist = np.linspace(-0.25*np.pi, 2.25*np.pi, 10000)
fig, ax = plt.subplots(1, 1, figsize=(3.5, 2.5))
ax.plot(philist, V(philist, lamb)/V(np.pi, lamb), label=r'$\lambda=$'+str(lamb), color='k')
ax.set_xlabel(r'$\varphi/\varphi_0$', labelpad = 5, fontsize=15)
ax.set_ylabel(r'$V(\varphi)/V(\varphi_{\rm fv})$', labelpad = 5, fontsize=15)
ax.xaxis.set_major_locator(plt.MultipleLocator(np.pi))
ax.xaxis.set_major_formatter(plt.FuncFormatter(multiple_formatter()))
a = ax.get_yticks().tolist()[1:-1:2]
ax.set_yticks(a)
a[0] = r'${:.0f}$'.format(a[0])
a[1] = r'${:.1f}$'.format(a[1])
a[2] = r'${:.0f}$'.format(a[2])
ax.set_yticklabels(a)
ax.grid(ls=':', color='darkgray', alpha=0.5)
# plt.legend()
plt.savefig('./plots/potential.pdf', dpi=500)
plt.tight_layout()
plt.show()
0.27428630597271203 0.03977359171321589 4.251836544352683 0.1 0.8999999999999999 4.047755453286628
[fluct_stdev(lambList[0], phi0List[0], tempList[ii]) for ii in [0,1,2,3]]
[0.27428630597271203, 0.30326206559910246, 0.33081957224035113, 0.3571216688274275]
[fluct_stdev_mom(lambList[0], phi0List[0], tempList[ii]) for ii in [0,1,2,3]]
[0.03977359171321589, 0.04522614119404094, 0.050656217225139674, 0.05606566994218313]
[phi0List[0]/fluct_stdev(lambList[0], phi0List[0], tempList[ii]) for ii in [0,1,2,3]]
[5.090532670392863, 4.6041478970906144, 4.220619088948682, 3.909769480468524]
[phi0List[0]/fluct_stdev_mom(lambList[0], phi0List[0], tempList[ii]) for ii in [0,1,2,3]]
[35.1052882440968, 30.872928017556298, 27.563514965790336, 24.904070584286945]
import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline
Vbare = lambda phi, lamb: 4.*nu* (- np.cos(phi) + 0.5 * lamb**2. * np.sin(phi)**2.)
Vdoub = lambda phi, lamb: 4.*nu* (np.cos(phi) + lamb**2. * np.cos(2.*phi) )
Vtrip = lambda phi, lamb: 4.*nu* (- np.sin(phi) - 4.*lamb**2. * np.cos(phi) * np.sin(phi) )
philist = np.linspace(0.5*np.pi, 1.5*np.pi, 10000)
yo = Vbare(philist, lambList[0])
xo = philist
plt.plot(philist, Vbare(philist, lambList[0]))
plt.plot(philist, Vdoub(philist, lambList[0]))
plt.plot(philist, Vtrip(philist, lambList[0]))
print(Vdoub(np.pi, lamb))
dyo = np.diff(yo,1)
dxo = np.diff(xo,1)
yfirst = dxo/dyo
xfirst = 0.5*(xo[:-1]+xo[1:])
dyfirst = np.diff(yfirst,1)
dxfirst = np.diff(xfirst,1)
ysecond = dyfirst/dxfirst
xsecond=0.5*(xfirst[:-1]+xfirst[1:])
#plt.plot(x[1:-1], xsecond)
0.01
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:14: RuntimeWarning: divide by zero encountered in true_divide
if True:
for tmp in [0,1,2,3]:
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
sigmamom = fluct_stdev_mom(lamb, phi0, temp)
print(sigmafld, sigmamom, right_Vmax.x, m2(lamb)**0.5, temp/m2(lamb)**0.5)
philist = np.linspace(0.75*np.pi, 2.25*np.pi, 10000)
fig, ax = plt.subplots(1, 1, figsize=(5, 2.5))
# ax.plot(philist, phi0**2.*V(philist, lamb), label=r'$\lambda=$'+str(lamb), color='k')
ax.plot(philist, V(philist, lamb)/V(np.pi, lamb), label=r'$\lambda=$'+str(lamb), color='k')
# ax.axvline(right_Vmax.x, ls='-', color='k', linewidth=0.5)
ax.axvline((4.*sigmafld+np.pi), ls='-', color='k', linewidth=0.5)
for nn in np.linspace(0, 5, 6):
ax.plot(np.pi + nn*sigmafld, V(np.pi + nn*sigmafld, lamb)/V(np.pi, lamb), \
marker='o', ms=3, color='r', label=int(nn))
ax.plot(np.pi + phi0/sigmafld, V(np.pi + phi0/sigmafld, lamb)/V(np.pi, lamb), \
marker='*', ms=3, color='g', label=int(nn))
ax.set_xlabel(r'$\varphi/\phi_0$', labelpad = 10)
ax.set_ylabel(r'$V(\varphi)/V(\varphi_{\rm fv})$', labelpad = 10)
ax.xaxis.set_major_locator(plt.MultipleLocator(np.pi))
ax.xaxis.set_major_formatter(plt.FuncFormatter(multiple_formatter()))
a = ax.get_xticks().tolist()[1:]
a[2] = a[1]
a[1] = (4.*sigmafld+np.pi)#right_Vmax.x
ax.set_xticks(a)
a[0] = r'$\pi$'
a[1] = r'${:.0f}$'.format((right_Vmax.x-np.pi)/sigmafld)+r'$\, \sigma_\varphi$'
a[2] = r'$2\pi$'
ax.set_xticklabels(a)
ax.grid(ls=':', color='darkgray', alpha=0.5)
plt.legend(bbox_to_anchor=(1., 1))
# plt.savefig('./plots/potential.pdf', dpi=500)
plt.tight_layout()
plt.show()
0.27428630597271203 0.03977359171321589 4.251836544352683 0.1 0.8999999999999999
0.30326206559910246 0.04522614119404094 4.251836544352683 0.1 1.0
0.33081957224035113 0.050656217225139674 4.251836544352683 0.1 1.0999999999999999
0.3571216688274275 0.05606566994218313 4.251836544352683 0.1 1.2
if False:
for tmp in range(len(tempList)):
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
def fun(y, x):
Φ, Π = y
#dydx = [Π, -Π/x+(np.sin(Φ) + 0.5*lamb**2.*np.sin(2.*Φ))*4.*nu]
# no friction term in thermal case
dydx = [Π, +(np.sin(Φ) + 0.5*lamb**2.*np.sin(2.*Φ))*4.*nu]
return dydx
#y0 = [5.69219082899233042025, 0.] # with friction
y0 = [4.823729806624866, 0.] # supercrit_instanton_file
# as far as fortran goes:
#4.8238 is super-critical
#4.823729806620144 is sub-critical
#y0 = [4.823729806625, 0.] # instanton_file
x = np.linspace(0, lenLat, 1000)
xplot = np.linspace(np.pi*0.8, 2.*np.pi, 1000)
sol = odeint(fun, y0, x)
plt.plot(x*np.sqrt(m2(lamb)), sol[:, 0], ls='-', label=r'$\phi(x)$')
plt.axhline(np.pi, ls=':', color='darkgray')
plt.legend(); plt.show()
plt.plot(xplot, Vinv(xplot, lamb), label='V', linewidth='1', alpha=0.3)
plt.plot(sol[:, 0], Vinv(sol[:, 0], lamb))
plt.legend(); plt.show()
ind = np.argmin(np.abs(sol[:,0] - np.pi)); print(ind, sol[ind,0])
inst = np.ones(nLat//2+1) * phieq
inst[:ind+1] = sol[:ind+1, 0]
inst = np.concatenate((inst[::-1],inst[1:-1]))
print(np.argmax(inst))
# np.save(supercrit_instanton_file(*exp_params), inst)
plt.plot(np.arange(nLat)-nLat//2, inst); plt.show()
print('Instanton saved.')
if True:
simLists = []
for tmp in [0,1,2,3]:
simLists.append([])
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
for sim in range(minSim, maxSim):
path2CLEANsim = clean_sim_location(*exp_params, sim)
path2RESTsim = rest_sim_location(*exp_params, sim)
if os.path.exists(path2RESTsim) and os.path.exists(path2CLEANsim):
simLists[-1].append(sim)
tp = 0# 0 for average, 1 for error
cp = 0
titl = [r'$\bar{\varphi}$', r'$\bar{\Pi}$', r'$\partial_x\bar{\varphi}$']
crit_rad = 40
crit_thresh = right_Vmax.x + 3.*sigmafld
win = 180
for tmp in [0,1,2,3]:
if tmp!=0: break
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
dats = random.sample(simLists[tmp], 6)
fig, ax = plt.subplots(2,3, figsize = (5.5,3.5), sharey=True, sharex=True, constrained_layout=True)
for si, sim in enumerate(dats):
ci,ti = np.divmod(si,3)
path2CLEANsim = clean_sim_location(*exp_params, sim)
fullreal, sim, tdecay, outcome = np.load(path2CLEANsim)
bubble2measure = fullreal[0]
nT, nN = np.shape(bubble2measure)
tcen, xcen = find_nucleation_center(bubble2measure, phieq, crit_thresh, crit_rad)
tcen -= 100
tl,tr = max(0, tcen-win), min(nT-1, tcen+win)
xl,xr = max(0, xcen-win), min(nN-1, xcen+win)
ext = np.array([xl-xcen,xr-xcen,tl-tcen-100,tr-tcen-100])*np.sqrt(m2(lamb))*dx
bubble2plot = bubble2measure[tl:tr,xl:xr]
im = ax[ci,ti].imshow(bubble2plot, interpolation=None, extent=ext, origin='lower', cmap='RdBu')
if False:
lavs = 6
nT, nN = np.shape(bubble2plot)
tt = np.linspace(tl-tcen, tr-tcen, nT)*np.sqrt(m2(lamb))*dx
xx = np.linspace(xl-xcen, xr-xcen, nN)*np.sqrt(m2(lamb))*dx
ttt1, xxx1 = np.meshgrid(tt, xx)
ax[ci,ti].contour(xxx1, ttt1, bubble2plot.T, levels=lavs, aspect='auto', interpolation=None, extent=ext, origin='lower', colors='k',linewidths=0.5)
if ci==1: ax[ci,ti].set(xlabel=r'$m \, r$')
if ti==0: ax[ci,ti].set(ylabel=r'$m \, t$')
ax[ci,ti].tick_params(direction='in', which='both', top=True, right=True)
ax[ci,ti].grid(ls=':', color='darkgray', alpha=0.5)
cbar = fig.colorbar(im, ax=ax[:, :], shrink=0.6, \
ticks=mticker.MultipleLocator(np.pi/2), \
format=mticker.FuncFormatter(multiple_formatter()))
cbar.ax.set_title(titl[0])
plt.savefig('./plots/examples_before'+str(temp)+'.pdf', bbox_inches='tight')
plt.show()
fig, ax = plt.subplots(2,3, figsize = (5.6,3.4), sharey=True, sharex=True, constrained_layout=True)
for si, sim in enumerate(dats):
ci,ti = np.divmod(si,3)
path2RESTsim = rest_sim_location(*exp_params, sim)
sim, fullreal, totbeta = np.load(path2RESTsim)
bubble2measure = fullreal[0]
nT, nN = np.shape(bubble2measure)
tcen, xcen = find_nucleation_center(bubble2measure, phieq, crit_thresh, crit_rad)
tcen -= 100
tl,tr = max(0, tcen-win), min(nT-1, tcen+win)
xl,xr = max(0, xcen-win), min(nN-1, xcen+win)
ext = np.array([xl-xcen,xr-xcen,tl-tcen-100,tr-tcen-100])*np.sqrt(m2(lamb))*dx
bubble2plot = bubble2measure[tl:tr,xl:xr]
im = ax[ci,ti].imshow(bubble2plot, interpolation=None, extent=ext, origin='lower', cmap='RdBu')
if False:
lavs = 6
nT, nN = np.shape(bubble2plot)
tt = np.linspace(tl-tcen, tr-tcen, nT)*np.sqrt(m2(lamb))*dx
xx = np.linspace(xl-xcen, xr-xcen, nN)*np.sqrt(m2(lamb))*dx
ttt1, xxx1 = np.meshgrid(tt, xx)
ax[ci,ti].contour(xxx1, ttt1, bubble2plot.T, levels=lavs, aspect='auto', interpolation=None, extent=ext, origin='lower', colors='k',linewidths=0.5)
if ci==1: ax[ci,ti].set(xlabel=r'$m \, r$')
if ti==0: ax[ci,ti].set(ylabel=r'$m \, t$')
ax[ci,ti].tick_params(direction='in', which='both', top=True, right=True)
ax[ci,ti].grid(ls=':', color='darkgray', alpha=0.5)
ax[ci,ti].legend(title=r'$v_{\rm COM}=$'+r'${:.2f}$'.format(totbeta), loc=4, fontsize=7, fancybox=True, frameon=True, framealpha=0.9, borderpad=0.15)
cbar = fig.colorbar(im, ax=ax[:, :], shrink=0.6, \
ticks=mticker.MultipleLocator(np.pi/2), \
format=mticker.FuncFormatter(multiple_formatter()))
cbar.ax.set_title(titl[0])
# plt.tight_layout()
plt.savefig('./plots/examples_after'+str(temp)+'.pdf', bbox_inches='tight')
plt.show()
No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument. No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument. No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument. No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument. No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument. No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
if True:
data_photo1 = np.load('./plots/bubble_data_original_2.npy')
rr1, ll1, rrwallfit1, llwallfit1, ttwallfit1, simulation1, ext1 = data_photo1
data_photo2 = np.load('./plots/bubble_data_deboosted_2.npy')
rr2, ll2, rrwallfit2, llwallfit2, ttwallfit2, simulation2, ext2 = data_photo2
data_photo3 = np.load('./plots/bubble_data_original_1.npy')
rr3, ll3, rrwallfit3, llwallfit3, ttwallfit3, simulation3, ext3 = data_photo3
data_photo4 = np.load('./plots/bubble_data_deboosted_1.npy')
rr4, ll4, rrwallfit4, llwallfit4, ttwallfit4, simulation4, ext4 = data_photo4
rrs = np.array([rr1, rr2, rr3, rr4])
lls = np.array([ll1, ll2, ll3, ll4])
rrwallfits = np.array([rrwallfit1, rrwallfit2, rrwallfit3, rrwallfit4])
llwallfits = np.array([llwallfit1, llwallfit2, llwallfit3, llwallfit4])
ttwallfits = np.array([ttwallfit1, ttwallfit2, ttwallfit3, ttwallfit4])
simulations = np.array([simulation1, simulation2, simulation3, simulation4])
exts = np.array([ext1, ext2, ext3, ext4])
exts = exts/dx2plot
exts[:,2:]+=20
exts = exts*dx2plot
ttwallfits+=20
dx2plot = dx * np.sqrt(m2(lamb))
fig, ax = plt.subplots(2, 2, figsize = (5.5, 5))
for ai, aa in enumerate(ax.flatten()):
rr, ll = rrs[ai], lls[ai]
rrwallfit, llwallfit, ttwallfit = rrwallfits[ai], llwallfits[ai], ttwallfits[ai]
simulation, ext = simulations[ai], exts[ai]
im0 = aa.imshow(simulation, interpolation='antialiased', extent=ext, origin='lower', cmap='RdBu', aspect='auto')
aa.plot(rr*dx2plot, ttwallfit*dx2plot, color='k', ls='--', linewidth=0.5)
aa.plot(ll*dx2plot, ttwallfit*dx2plot, color='k', ls='--', linewidth=0.5)
aa.plot(rrwallfit*dx2plot, ttwallfit*dx2plot, color='k', ls='-', linewidth=0.7)
aa.plot(llwallfit*dx2plot, ttwallfit*dx2plot, color='k', ls='-', linewidth=0.7)
cbar = plt.colorbar(im0, ax=aa, shrink=0.9, \
ticks=mticker.MultipleLocator(np.pi/2), \
format=mticker.FuncFormatter(multiple_formatter()))
cbar.ax.set_title(r'$\bar{\varphi}$')
aa.tick_params(direction='in', which='both', top=True, right=True)
aa.grid(ls=':', color='darkgray', alpha=0.5)
aa.set_ylabel(r'$m \, t$', labelpad = 1)
aa.set_xlabel(r'$m \, r$', labelpad = 1)
aa.set_aspect(1.5)
plt.tight_layout()
plt.savefig('./plots/examples_deboost.pdf', dpi=500)
plt.show()
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:11: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray. # This is added back by InteractiveShellApp.init_path() /cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:12: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray. if sys.path[0] == "": /cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:13: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray. del sys.path[0] /cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:14: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray. /cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:15: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray. from ipykernel import kernelapp as app
tp = 0# 0 for average, 1 for error
titl = [r'$\bar{\varphi}$', r'$\bar{\Pi}$', r'$\partial_r\bar{\varphi}$']
for ii, average_bubble in enumerate(bubbleList):
fig, ax = plt.subplots(1, 3, figsize = (8,2.5))
for cp in range(3): # 0 - field, 1 - momentum, 2 - gradient
tmp = tmpList[ii]
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
crit_rad = 20
crit_thresh = right_Vmax.x + 2.*sigmafld
win = 120
bubble2measure = average_bubble[0,0]
nT, nN = np.shape(bubble2measure)
tcen, xcen = find_nucleation_center(bubble2measure, phieq, crit_thresh, crit_rad)
xcen+=1
tl,tr = max(0, tcen-win), min(nT-1, tcen+win)
xl,xr = max(0, xcen-win), min(nN-1, xcen+win)
ext = np.array([xl-xcen,xr-xcen,tl-tcen,tr-tcen])*np.sqrt(m2(lamb))*dx
bubble2plot = average_bubble[tp,cp][tl:tr,xl:xr]
im = ax[cp].imshow(bubble2plot, interpolation=None, extent=ext, origin='lower', cmap='RdBu')
if cp==0: cbar = plt.colorbar(im, ax=ax[cp], shrink=0.8, \
ticks=mticker.MultipleLocator(np.pi/2), \
format=mticker.FuncFormatter(multiple_formatter()))
else:
cbar = plt.colorbar(im, ax=ax[cp], shrink=0.8)
lls = cbar.ax.get_ylim()
a = cbar.ax.get_yticks().tolist()
a = [round(al,2) for ja,al in enumerate(a) if ja%2==1][::]
cbar.ax.set_yticks(a)
#cbar.ax.set_ylim(a[0], a[-1])
cbar.ax.set_ylim(lls)
a = [r'${:.1f}$'.format(al) for al in a]
cbar.ax.set_yticklabels(a)
cbar.ax.set_title(titl[cp])
nT, nN = np.shape(bubble2plot)
tt = np.linspace(tl-tcen, tr-tcen, nT)*np.sqrt(m2(lamb))*dx
xx = np.linspace(xl-xcen, xr-xcen, nN)*np.sqrt(m2(lamb))*dx
ttt1, xxx1 = np.meshgrid(tt, xx)
lavs = [8, 6, 8][cp]
ax[cp].grid(True, color='k', linewidth=0.5, ls=':', alpha=0.2)
ax[cp].contour(xxx1, ttt1, bubble2plot.T, levels=lavs, aspect='auto', interpolation=None, extent=ext, origin='lower', colors='k',linewidths=0.5)
labss = r'$T/m = {}$'.format(round(temp/np.sqrt(m2(lamb)),1))
ax[cp].grid(ls=':', color='darkgray', alpha=0.7)
ax[cp].xaxis.set_minor_locator(MultipleLocator(1))
ax[cp].yaxis.set_minor_locator(MultipleLocator(1))
ax[cp].set_xlabel(r'$m \, r$', labelpad = 1)
ax[0].set_ylabel(r'$m \, t$', labelpad = 1)
ax[cp].yaxis.set_ticks_position('both')
ax[cp].xaxis.set_ticks_position('both')
ax[cp].tick_params(which='both', axis="y", direction="in")
ax[cp].tick_params(which='both', axis="x", direction="in")
plt.tight_layout()
plt.savefig('./plots/average_bubble'+str(temp)+'.pdf', bbox_inches='tight')
plt.show()
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:50: UserWarning: The following kwargs were not used by contour: 'aspect', 'interpolation'
if False:
for ii, average_bubble in enumerate(bubbleList):
average_bubble = bubbleList[ii]
tmp = tmpList[ii]
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
print(*exp_params)
critical_sim = extract_data(nLat, critical_sim_file(*exp_params))
precursor_sim = extract_data(nLat, precursor_sim_file(*exp_params))
titls = [r'$\rm Critical$', r'$\rm Precursor$']
tdecay = max(0, np.shape(critical_sim)[1] - nLat//2)
crit_thresh = right_Vmax.x + 5.*sigmafld
critical_sim, crit_rad = centre_bubble(critical_sim, tdecay, phieq, crit_thresh)
critical_sim = remove_collisions(critical_sim, phieq)
print('Decay time standard, radius:', tdecay, crit_rad)
critical_sim = critical_sim[0]
precursor_sim = precursor_sim[0]
# critical_sim = gaussian_filter(critical_sim, sigma=5, mode='nearest')
# precursor_sim = gaussian_filter(precursor_sim, sigma=5, mode='nearest')
truebubble = critical_sim[:nLat//4, 3*nLat//8 : 5*nLat//8]
precursor = precursor_sim[:nLat//4, 3*nLat//8 : 5*nLat//8]
print(np.shape(precursor))
nN, nT = np.shape(truebubble)
tcen, xcen = find_nucleation_center(truebubble, phieq, crit_thresh, crit_rad//3)
xcen += 1
tl,tr,xl,xr = -tcen, nT-1-tcen, -xcen, nN-1-xcen
exts1 = np.array([xl, xr, tl, tr])*dx*np.sqrt(m2(lamb))
tt1, xx1 = np.linspace(tl, tr, nT)*dx*np.sqrt(m2(lamb)), np.linspace(xl, xr, nN)*dx*np.sqrt(m2(lamb))
ttt1, xxx1 = np.meshgrid(tt1, xx1)
tl,tr,xl,xr = 0, nT-1, -xcen, nN-1-xcen
exts2 = np.array([xl, xr, tl, tr])*dx*np.sqrt(m2(lamb))
tt2, xx2 = np.linspace(tl, tr, nT)*dx*np.sqrt(m2(lamb)), np.linspace(xl, xr, nN)*dx*np.sqrt(m2(lamb))
ttt2, xxx2 = np.meshgrid(tt2, xx2)
fig, ax = plt.subplots(1, 2, figsize = (6,3))
ax[0].contour(xxx1, ttt1, truebubble.T, levels=3, aspect='auto', interpolation=None, extent=exts1, origin='lower', colors='k', linewidths=0.5)
im0 = ax[0].imshow(truebubble, aspect='auto', interpolation='none', extent=exts1, origin='lower', cmap='RdBu')
clb0 = plt.colorbar(im0, ax = ax[0], ticks = mticker.MultipleLocator(np.pi/2), format = mticker.FuncFormatter(multiple_formatter()))
clb0.ax.set_title(r'$\bar{\varphi}$')
ax[1].contour(xxx2, ttt2, np.abs(np.pi-precursor.T), levels=3, aspect='auto', interpolation=None, extent=exts2, origin='lower', colors='k', linewidths=0.5)
im1 = ax[1].imshow(precursor, aspect='auto', interpolation='none', extent=exts2, origin='lower', cmap='RdBu')
clb1 = plt.colorbar(im1, ax = ax[1], ticks = mticker.MultipleLocator(np.pi/2), format = mticker.FuncFormatter(multiple_formatter()))
clb1.ax.set_title(r'$\bar{\varphi}$')
for ai, aa in enumerate(ax):
# aa.text(5.5,.5,'TEST TEST TEST TEST',
# bbox={'facecolor':'white','alpha':1,'edgecolor':'none','pad':1},
# ha='center', va='center')
aa.text([6.8,6.2][ai], [12.,18.][ai], titls[ai], ha='center', va='center', \
bbox={'boxstyle':'round','facecolor':'white','alpha':0.85,'edgecolor':'k','pad':0.3}, fontsize=10)
aa.grid(ls=':', color='darkgray', alpha=0.7)
aa.xaxis.set_minor_locator(MultipleLocator(1))
aa.yaxis.set_minor_locator(MultipleLocator(1))
# aa.set_xlabel(r'$\phi_0^{-1} \sqrt{V_0} \; r$')
# aa.set_ylabel(r'$\phi_0^{-1} \sqrt{V_0} \; t$')
aa.set_xlabel(r'$m \, r$')
aa.set_ylabel(r'$m \, t$')
aa.yaxis.set_ticks_position('both')
aa.xaxis.set_ticks_position('both')
aa.tick_params(which='both', axis="y", direction="in")
aa.tick_params(which='both', axis="x", direction="in")
# a = aa.get_xticks().tolist()[1:-1:]
# a = [round(al,2) for al in a]
# aa.set_xticks(a)
# aa.set_xticklabels(a)
a = aa.get_yticks().tolist()[1:-1:2]
a = [round(al,2) for al in a]
aa.set_yticks(a)
a = [r'${:.0f}$'.format(al) for al in a]
aa.set_yticklabels(a)
plt.tight_layout()
plt.savefig('./plots/critical_and_precursor'+str(temp)+'.pdf', bbox_inches='tight', dpi=500)
plt.show()
if False:
titls = [r'$\rm Critical$', r'$\rm Subcritical$']
super_instanton_sim = extract_data(nLat, supercrit_instanton_file(*exp_params))
instanton_sim = extract_data(nLat, instanton_file(*exp_params))
tdecay = max(0, np.shape(super_instanton_sim)[1] - nLat//2)
crit_thresh = right_Vmax.x + 5.*sigmafld
super_instanton_sim, crit_rad = centre_bubble(super_instanton_sim, tdecay, phieq, crit_thresh)
super_instanton_sim = remove_collisions(super_instanton_sim, phieq)
for tmp in range(len(tempList)):
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
# np.save(supercrit_instanton_file(*exp_params), super_instanton_sim[0,0,:])
# np.save(instanton_file(*exp_params), instanton_sim[0,0,:])
print('Decay time standard, radius:', tdecay, crit_rad)
super_bubble = super_instanton_sim[0, 40:nLat//3+41, nLat//3:2*nLat//3+1]
bubble = instanton_sim[0, 40:nLat//3+41, nLat//3:2*nLat//3+1]
nN, nT = np.shape(super_bubble)
tcen, xcen = find_nucleation_center(super_bubble, phieq, crit_thresh, crit_rad//3)
tl,tr,xl,xr = -tcen, nT-1-tcen, -xcen, nN-1-xcen
exts1 = np.array([xl, xr, tl, tr])*dx*np.sqrt(m2(lamb))
tt1, xx1 = np.linspace(tl, tr, nT)*dx*np.sqrt(m2(lamb)), np.linspace(xl, xr, nN)*dx*np.sqrt(m2(lamb))
ttt1, xxx1 = np.meshgrid(tt1, xx1)
tl,tr,xl,xr = 0, nT-1, -xcen, nN-1-xcen
exts2 = np.array([xl, xr, tl, tr])*dx*np.sqrt(m2(lamb))
tt2, xx2 = np.linspace(tl, tr, nT)*dx*np.sqrt(m2(lamb)), np.linspace(xl, xr, nN)*dx*np.sqrt(m2(lamb))
ttt2, xxx2 = np.meshgrid(tt2, xx2)
fig, ax = plt.subplots(1, 2, figsize = (6,3))
ax[0].contour(xxx1, ttt1, super_bubble.T, levels=4, aspect='auto', interpolation='none', extent=exts1, origin='lower', colors='k', linewidths=0.5)
im0 = ax[0].imshow(super_bubble, aspect='auto', interpolation='none', extent=exts1, origin='lower', cmap='RdBu')
clb0 = plt.colorbar(im0, ax = ax[0], ticks = mticker.MultipleLocator(np.pi/2), format = mticker.FuncFormatter(multiple_formatter()))
clb0.ax.set_title(r'$\bar{\varphi}$')
ax[1].contour(xxx2, ttt2, bubble.T, levels=7, aspect='auto', interpolation='none', extent=exts2, origin='lower', colors='k', linewidths=0.5)
im1 = ax[1].imshow(bubble, aspect='auto', interpolation='none', extent=exts2, origin='lower', cmap='RdBu')
clb1 = plt.colorbar(im1, ax = ax[1], ticks = mticker.MultipleLocator(np.pi/2), format = mticker.FuncFormatter(multiple_formatter()))
clb1.ax.set_title(r'$\bar{\varphi}$')
for ai, aa in enumerate(ax):
# aa.text(5.5,.5,'TEST TEST TEST TEST',
# bbox={'facecolor':'white','alpha':1,'edgecolor':'none','pad':1},
# ha='center', va='center')
aa.text([9.,7.2][ai], [-7.8,2.][ai], titls[ai], ha='center', va='center', \
bbox={'boxstyle':'round','facecolor':'white','alpha':0.85,'edgecolor':'k','pad':0.3}, fontsize=9)
aa.grid(ls=':', color='darkgray', alpha=0.7)
aa.xaxis.set_minor_locator(MultipleLocator(1))
aa.yaxis.set_minor_locator(MultipleLocator(1))
# aa.set_xlabel(r'$\phi_0^{-1} \sqrt{V_0} \; r$')
# aa.set_ylabel(r'$\phi_0^{-1} \sqrt{V_0} \; t$')
aa.set_xlabel(r'$m \, r$')
aa.set_ylabel(r'$m \, t$')
aa.yaxis.set_ticks_position('both')
aa.xaxis.set_ticks_position('both')
aa.tick_params(which='both', axis="y", direction="in")
aa.tick_params(which='both', axis="x", direction="in")
plt.tight_layout()
plt.savefig('./plots/bare_spheleron.pdf', dpi=500)
plt.show()
tmp = 0
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
print(*exp_params)
instanton = np.load(instanton_file(*exp_params))# - np.pi
supercritinstanton = np.load(supercrit_instanton_file(*exp_params))# - np.pi
print(np.shape(instanton), np.shape(supercritinstanton))
xmin, xmax = nLat//3, nLat*2//3
#plt.plot(xlist[xmin:xmax], instanton[xmin:xmax])
#plt.plot(xlist[xmin:xmax], supercritinstanton[xmin:xmax])
plt.plot(xlist[xmin:xmax], supercritinstanton[xmin+1:xmax+1] - instanton[xmin:xmax])
plt.plot(xlist[xmin:xmax], supercritinstanton[xmin:xmax] - super_instanton_sim[0,0,xmin:xmax])
plt.plot(xlist[xmin:xmax], instanton[xmin:xmax] - instanton_sim[0,0,xmin:xmax])
plt.show()
crit_soln_f90 = '(/'
for ind, iii in enumerate(instanton):
crit_soln_f90 = crit_soln_f90 + str(iii)
if ind != len(instanton)-1:
crit_soln_f90 = crit_soln_f90 + ', '
crit_soln_f90 += '/)'
#print(crit_soln_f90)
if False:
tmp = 5
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
undecayed_sims = np.load(sims_notdecayed_file(*exp_params, minSim, maxSim, nTimeMAX))
print(undecayed_sims[:,0])
decay_times = np.load(decay_times_file(*exp_params, minSim, maxSim, nTimeMAX))
minDecTime = 4*nLat
alltimes = decay_times[:,1]
simList2Do = decay_times[alltimes>=minDecTime, 0]
print(simList2Do)
# sims_decayed = np.load(sims_decayed_file(*exp_params, minSim, maxSim, nTimeMAX))
choose = random.sample(undecayed_sims[:,0].tolist(), 1)
for sim in choose:
path2sim = sim_location(*exp_params, sim)
real, cho = get_realisation(nLat, sim, phieq, path2sim)
nC, nT, nN = np.shape(real)
print('Sim, nT, outcome', sim, nT, cho)
bubble = real[0,:]#nT-3*nLat//4]
nT, nN = np.shape(bubble)
ax1 = simple_imshow(bubble, [0,nT,0,nN], 'Simulation '+str(sim))
plt.show()
fftbubble = np.fft.fft2(bubble)
nT, nN = np.shape(fftbubble)
#fftbubble = fftbubble[:nT//2,:nN//2]
fftbubble[:, nN*4//10:nN*6//10+1] = 0.
ax2 = simple_imshow(np.log(fftbubble.real), [0,nT,0,nN], 'Simulation '+str(sim))
bubble2 = np.fft.ifft2(fftbubble)
nT, nN = np.shape(bubble2)
ax3 = simple_imshow(bubble2.real, [0,nT,0,nN], 'Simulation '+str(sim))
ax4 = simple_imshow(bubble2.imag, [0,nT,0,nN], 'Simulation '+str(sim))
ax5 = simple_imshow(bubble2.real-bubble, [0,nT,0,nN], 'Simulation '+str(sim))
if False:
# if sim in sims_decayed[:,0]:
path2CLEANsim = clean_sim_location(*exp_params, sim)
real, sim, tdecay, outcome = np.load(path2CLEANsim)
bubble = real[0,-nLat:]
nT, nN = np.shape(bubble)
simple_imshow(bubble, [0,nT,0,nN], 'Simulation '+str(sim))
print('Sim, outcome0, outcome1, tdecay', sim, outcome, cho, tdecay)
# Classify decays
fig, ax = plt.subplots(1, len(tempList), figsize = (15,2.5))
plt.style.use('seaborn-whitegrid') # nice and clean grid
cols = cycle(['purple', 'forestgreen', 'orange', 'blue'])
for tmp in range(len(tempList)):
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
undecayed_sims = np.load(sims_notdecayed_file(*exp_params, minSim, maxSim, nTimeMAX))
decayed_sims = np.load(sims_decayed_file(*exp_params, minSim, maxSim, nTimeMAX))
decay_times = np.load(decay_times_file(*exp_params, minSim, maxSim, nTimeMAX))
outcomes = decayed_sims[:,1]; print(len(outcomes))
labs = r'$T/m={:.1f}$'.format(temp/np.sqrt(m2(lamb)))
ax[tmp].hist(outcomes, density=True, bins=2, rwidth=0.8, color=next(cols), alpha=0.9, linewidth=0.5)
ax[tmp].set_xlabel(r'$2\pi \quad\quad 0$')
ax[tmp].set_ylabel('PDF')
ax[tmp].set_title(labs)
plt.tight_layout()
plt.savefig('./plots/vacuum_choice.pdf', rasterize=True); plt.show()
1003 2192 3165 3711 3910 3984
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:22: MatplotlibDeprecationWarning: savefig() got unexpected keyword argument "rasterize" which is no longer supported as of 3.3 and will become an error in 3.6
tmp = tmpList[0]
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
instanton = instantonList[tmp]
gradinst = (np.roll(instanton,-1) - instanton)/dx2plot
testgrd = fldcritList[tmp][-1,:]*dx/dx2plot
testfld = fldcritList[tmp][0,:]
testfld = (np.roll(testfld,-1) - testfld)/dx2plot
xtest = (np.arange(nLat) - nLat//2) * dx2plot
cds = np.array(np.linspace(nLat//2-len(testgrd)//2, nLat//2+len(testgrd)//2, len(testgrd)), dtype='int')
fig, ax = plt.subplots(1, 1, figsize = (5,3.5))
ax.plot(xtest[cds], testgrd, ls='-', label='gradient of average bubble')
ax.plot(xtest[cds], testfld, ls='-', label='average gradient of average bubble')
ax.plot(xtest[cds], gradinst[cds], ls='-',color='k', label='gradient of critical lattice solution')
ax.legend()
ax.set_ylabel(r'$\partial_x \phi_{\rm crit}(x)$')
ax.set_xlabel(r'$\phi_0^{-1} \sqrt{V_0} \; x$')
plt.grid(True, ls=':', color='lightgray', alpha=0.5)
plt.show()
if True:
for tmp in range(len(tempList)):
if tmp==4: break
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
# 0 for field; 1 for momentum
find = 0
# modes to plot:
aa, bb = 1, knyq-1
# time steps to plot
tstep = 10
# if tind not in np.array([0,2,7,27,46]): continue
ALL_powspec1 = np.load(powspec_tlist_file(*exp_params, minSim, 2000))
tlist, PSfld1 = ALL_powspec1[0][::tstep], ALL_powspec1[1][:, find, ::tstep, aa:bb]
del ALL_powspec1
ALL_powspec2 = np.load(powspec_tlist_file(*exp_params, 2000, maxSim))
tlist, PSfld2 = ALL_powspec2[0][::tstep], ALL_powspec2[1][:, find, ::tstep, aa:bb]
del ALL_powspec2
fig, ax = plt.subplots(1,1, figsize = (6,4.))
plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), pspec(lamb, phi0, temp)[aa:bb], color='k', linewidth=0.5)
for tind, tt in enumerate(tlist):
curve = np.nanmean(np.concatenate((PSfld1[:,tind], PSfld2[:,tind]), axis=0), axis=0)
lab = r'${:.1f}$'.format(round(tt*dx*np.sqrt(m2(lamb)),1))
plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), curve, ls='-', label=lab)
del PSfld1, PSfld2
h, l = ax.get_legend_handles_labels() # Extracting handles and labels
handles = [plt.plot([],marker="", ls="")[0]] + h
labels = [r'$m\,t:$'] + l # Merging labels
leg = ax.legend(handles, labels, ncol=2, frameon=False, loc=1)
ax.grid(which='both', color='darkgray', ls=':', alpha=0.4)
ax.tick_params(direction='in', which='both', top=True, right=True)
ax.set_xscale('log')
ax.set_ylabel(r'$\left\langle \left|\bar{\varphi}_k(t)\right|^2 \right\rangle $')
ax.set_xlabel(r'$k/m$')
plt.tight_layout()
plt.savefig('./plots/powespec_tevol.pdf')
plt.show()
if True:
for tmp in range(len(tempList)):
if tmp==4: break
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
# 0 for field; 1 for momentum
find = 1
# modes to plot:
aa, bb = 1, knyq-1
# time steps to plot
tstep = 10
ALL_powspec1 = np.load(powspec_tlist_file(*exp_params, minSim, 2000))
tlist, PSfld1 = ALL_powspec1[0][::tstep], ALL_powspec1[1][:, find, ::tstep, aa:bb]
del ALL_powspec1
ALL_powspec2 = np.load(powspec_tlist_file(*exp_params, 2000, maxSim))
tlist, PSfld2 = ALL_powspec2[0][::tstep], ALL_powspec2[1][:, find, ::tstep, aa:bb]
del ALL_powspec2
fig, ax = plt.subplots(1,1, figsize = (4.7,3.))
plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), (w2(lamb)*pspec(lamb, phi0, temp))[aa:bb], color='k', linewidth=0.5)
for tind, tt in enumerate(tlist):
curve = np.nanmean(np.concatenate((PSfld1[:,tind], PSfld2[:,tind]), axis=0), axis=0)
lab = r'${:.1f}$'.format(round(tt*dx*np.sqrt(m2(lamb)),1))
plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), curve, ls='-', label=lab)
del PSfld1, PSfld2
h, l = ax.get_legend_handles_labels() # Extracting handles and labels
handles = [plt.plot([],marker="", ls="")[0]] + h
labels = [r'$m\,t:$'] + l # Merging labels
leg = ax.legend(handles, labels, ncol=3, frameon=False, bbox_to_anchor=(1.5, 0.5))
ax.grid(which='both', color='darkgray', ls=':', alpha=0.4)
ax.tick_params(direction='in', which='both', top=True, right=True)
ax.set_xscale('log')
ax.set_ylabel(r'$\left\langle \left|\bar{\Pi}_k(t)\right|^2 \right\rangle $')
ax.set_xlabel(r'$k/m$')
plt.tight_layout()
plt.savefig('./plots/powespec_tevol.pdf')
plt.show()
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:43: UserWarning: Tight layout not applied. The left and right margins cannot be made large enough to accommodate all axes decorations.
if True:
savemeff = np.zeros((4))
for tmp in range(len(tempList)):
if tmp==4: break
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
# modes to plot:
aa, bb = 1, knyq//5
ALL_powspec1 = np.load(powspec_tlist_file(*exp_params, minSim, 2000))
tlist, PSfld1, PSmom1 = ALL_powspec1[0], ALL_powspec1[1][:, 0, -50:, aa:bb], ALL_powspec1[1][:, 1, -50:, aa:bb]
del ALL_powspec1
ALL_powspec2 = np.load(powspec_tlist_file(*exp_params, 2000, maxSim))
tlist, PSfld2, PSmom2 = ALL_powspec2[0], ALL_powspec2[1][:, 0, -50:, aa:bb], ALL_powspec2[1][:, 1, -50:, aa:bb]
del ALL_powspec2
avPSfld = np.nanmean(np.nanmean(np.concatenate((PSfld1, PSfld2), axis=0), axis=1), axis=0)
avPSmom = np.nanmean(np.nanmean(np.concatenate((PSmom1, PSmom2), axis=0), axis=1), axis=0)
del PSfld1, PSfld2, PSmom1, PSmom2
curve = (avPSmom/avPSfld)
norml = 1./ phi0 / np.sqrt(2. * lenLat)
w2m = lambda ks, m: ks**2. + m**2.
pofk = lambda ks, m: w2m(ks, m)
pred_fit = lambda x, data: sco.curve_fit(pofk, x, data, p0=0.09)[0]
fig, ax = plt.subplots(1,1, figsize = (4.7,3.))
plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), pofk(klist[aa:bb], np.sqrt(m2(lamb))), color='k', lw=0.5, ls='-')
lab1 = r'$m={:.4f}$'.format(np.sqrt(m2(lamb)), temp)
plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), curve, ls='-', label=lab1)
best_ps = pred_fit(klist[aa:bb], curve)
lab2 = r'$m={:.4f}$'.format(*best_ps)
savemeff[tmp] = best_ps[0]
plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), pofk(klist[aa:bb], *best_ps), color='r', lw=1, ls=':', label=lab2)
leg = ax.legend(ncol=1, frameon=False, loc='best')#, bbox_to_anchor=(2.3, 1.))
ax.grid(which='both', color='darkgray', ls=':', alpha=0.4)
ax.tick_params(direction='in', which='both', top=True, right=True)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_ylabel(r'$\left\langle \left|\bar{\Pi}_k(t)\right|^2/\left|\bar{\varphi}_k(t)\right|^2 \right\rangle $')
ax.set_xlabel(r'$k/m$')
plt.tight_layout()
#plt.savefig('./plots/bestfit_masses.pdf')
plt.show()
print(savemeff)
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:19: RuntimeWarning: Mean of empty slice /cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:20: RuntimeWarning: Mean of empty slice
[0.08631948 0.08449074 0.08163867 0.08173309]
if True:
norml = 1./ phi0 / np.sqrt(2. * lenLat)
w2m = lambda ks, msq: ks**2. + msq
pofk0 = lambda ks, msq, te: norml / w2m(ks, msq)**0.25 * np.sqrt(2./(np.exp(w2m(ks, msq)**0.5/te) - 1.))
#pofk = lambda ks, msq, te: norml / w2m(ks, msq)**0.25 * np.sqrt(2./(w2m(ks, msq)**0.5/te + 0.5*(w2m(ks, msq)**0.5/te)**2. ))
pofk = lambda ks, msq, te: norml / w2m(ks, msq)**0.25 * np.sqrt(2./(w2m(ks, msq)**0.5/te))
f_pred = lambda ks, msq, te: np.abs(pofk(ks, msq, te))**2.
f_pred1 = lambda ks, te: np.abs(pofk(ks, savemeff[tmp]**2., te))**2.
f_pred2 = lambda ks, te: np.abs(pofk0(ks, m2(lamb), te))**2.
pred_fit = lambda x, data: sco.curve_fit(f_pred1, x, data)[0]
save_teff_from_fld = np.zeros((4))
for tmp in range(len(tempList)):
if tmp==4: break
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
# modes to plot:
aa, bb = 1, knyq//5
# time slices to average
am, bm = -50, -1
ALL_powspec1 = np.load(powspec_tlist_file(*exp_params, minSim, 2000))
tlist, PSfld1 = ALL_powspec1[0], ALL_powspec1[1][:, 0, am:bm, aa:bb]
del ALL_powspec1
ALL_powspec2 = np.load(powspec_tlist_file(*exp_params, 2000, maxSim))
tlist, PSfld2 = ALL_powspec2[0], ALL_powspec2[1][:, 0, am:bm, aa:bb]
del ALL_powspec2
PSfld = np.concatenate((PSfld1, PSfld2), axis=0)
curve = np.nanmean(np.nanmean(PSfld, axis=1), axis=0)
del PSfld1, PSfld2
fig, ax = plt.subplots(1,1, figsize = (4.7,3.))
plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), f_pred2(klist[aa:bb], temp), color='k', lw=1, ls='-', label=r'$\rm BE, \; initially$')
lab1 = r'${:.1f} < m\,t $'.format(round(tlist[am]*dx*np.sqrt(m2(lamb)),1)) + r'$< {:.1f}$'.format(round(tlist[bm]*dx*np.sqrt(m2(lamb)),1))
[plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), ii, ls='-', label=(lab1 if i0==0 else None)) for i0, ii in enumerate(np.nanmean(PSfld, axis=0))]
#plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), curve, ls='-', lw=2, label=lab1)
best_ps = pred_fit(klist[aa:bb], curve)
lab2 = r'$T={:.4f}$'.format(*best_ps)
save_teff_from_fld[tmp] = best_ps[-1]
plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), f_pred1(klist[aa:bb], *best_ps), color='r', lw=1, ls=':', label=lab2)
leg = ax.legend(title=r'$T={:.2f}$'.format(temp), ncol=1, frameon=False, loc='best')
ax.grid(which='both', color='darkgray', ls=':', alpha=0.4)
ax.tick_params(direction='in', which='both', top=True, right=True)
ax.set_xscale('log')
ax.set_ylabel(r'$\left\langle \left|\bar{\varphi}_k(t)\right|^2 \right\rangle $')
ax.set_xlabel(r'$k/m$')
plt.tight_layout()
#plt.savefig('./plots/powespec_tevol.pdf')
plt.show()
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:32: RuntimeWarning: Mean of empty slice
if True:
norml = 1./ phi0 / np.sqrt(2. * lenLat)
w2m = lambda ks, msq: ks**2. + msq
pofk0 = lambda ks, msq, te: norml * w2m(ks, msq)**0.25 * np.sqrt(2./(np.exp(w2m(ks, msq)**0.5/te) - 1.))
#pofk = lambda ks, msq, te: norml * w2m(ks, msq)**0.25 * np.sqrt(2./(w2m(ks, msq)**0.5/te + 0.5*(w2m(ks, msq)**0.5/te)**2. ))
pofk = lambda ks, msq, te: norml * w2m(ks, msq)**0.25 * np.sqrt(2./(w2m(ks, msq)**0.5/te))
f_pred = lambda ks, msq, te: np.abs(pofk(ks, msq, te))**2.
f_pred1 = lambda ks, te: np.abs(pofk(ks, savemeff[tmp]**2., te))**2.
f_pred2 = lambda ks, te: np.abs(pofk0(ks, m2(lamb), te))**2.
pred_fit = lambda x, data: sco.curve_fit(f_pred1, x, data)[0]
save_teff_from_mom = np.zeros((4))
for tmp in range(len(tempList)):
if tmp==4: break
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
# modes to plot:
aa, bb = 1, knyq//5
# time slices to average
am, bm = -50, -1
ALL_powspec1 = np.load(powspec_tlist_file(*exp_params, minSim, 2000))
tlist, PSfld1 = ALL_powspec1[0], ALL_powspec1[1][:, 1, am:bm, aa:bb]
del ALL_powspec1
ALL_powspec2 = np.load(powspec_tlist_file(*exp_params, 2000, maxSim))
tlist, PSfld2 = ALL_powspec2[0], ALL_powspec2[1][:, 1, am:bm, aa:bb]
del ALL_powspec2
PSfld = np.concatenate((PSfld1, PSfld2), axis=0)
curve = np.nanmean(np.nanmean(PSfld, axis=1), axis=0)
del PSfld1, PSfld2
fig, ax = plt.subplots(1,1, figsize = (4.7,3.))
#plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), f_pred2(klist[aa:bb], temp), color='k', lw=1, ls='-', label=r'$\rm BE, \; initially$')
lab1 = r'${:.1f} < m\,t $'.format(round(tlist[am]*dx*np.sqrt(m2(lamb)),1)) + r'$< {:.1f}$'.format(round(tlist[bm]*dx*np.sqrt(m2(lamb)),1))
[plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), ii, ls='-', label=(lab1 if i0==0 else None)) for i0, ii in enumerate(np.nanmean(PSfld, axis=0))]
#plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), curve, ls='-', lw=2, label=lab1)
best_ps = pred_fit(klist[aa:bb//10], curve[:bb//10-1])
lab2 = r'$T={:.4f}$'.format(*best_ps)
save_teff_from_mom[tmp] = best_ps[-1]
plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), f_pred1(klist[aa:bb], *best_ps), color='r', lw=1, ls=':', label=lab2)
lab3 = r'$T={:.4f}$'.format(save_teff_from_fld[tmp])
plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), f_pred1(klist[aa:bb], save_teff_from_fld[tmp]), color='g', lw=1, ls=':', label=lab3)
leg = ax.legend(title=r'$T={:.2f}$'.format(temp), ncol=1, frameon=False, loc='best')
ax.grid(which='both', color='darkgray', ls=':', alpha=0.4)
ax.tick_params(direction='in', which='both', top=True, right=True)
ax.set_xscale('log')
ax.set_ylabel(r'$\left\langle \left|\bar{\Pi}_k(t)\right|^2 \right\rangle $')
ax.set_xlabel(r'$k/m$')
plt.tight_layout()
#plt.savefig('./plots/powespec_tevol.pdf')
plt.show()
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:32: RuntimeWarning: Mean of empty slice
fig, ax = plt.subplots(1,1, figsize = (4.7,3.))
plt.plot(tempList[tmpList], save_teff_from_fld, label=r'$\rm Best \; Fit \; RJ \; from Field$')
plt.plot(tempList[tmpList], save_teff_from_mom, label=r'$\rm Best \; Fit \; RJ \; from Momentum$')
plt.plot(tempList[tmpList], tempList[tmpList], label=r'$\rm Bare \; Temperature$')
plt.legend()
plt.show()
exp_bounce = lambda x, a, b: a * np.sqrt(b/x) * np.exp(-b/x)
get_best_bounce = lambda x, y: sco.curve_fit(exp_bounce, x, y)[0]
fig, ax = plt.subplots(1, 3, figsize = (9.7,3.), gridspec_kw={'width_ratios': [4, 2.3, 2.3]})
gammas1 = np.zeros((len(tempList), 2))
cls = cycle(allcolors[:4])
for tmp in range(len(tempList[:4])):
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
labss = r'${}$'.format(round(temp/np.sqrt(m2(lamb)),1))
lmmax, lmmin = 5*nLat, 1.5*nLat#2*nLat/3
decay_times = np.load(decay_times_file(*exp_params, minSim, maxSim, nTimeMAX))
ndcys = maxSim - minSim
decay_times = np.sort(decay_times[:,1])
decay_times = decay_times[decay_times < lmmax]
ax[0].fill_betweenx(np.linspace(-1,2,10), lmmin*dx*np.sqrt(m2(lamb)), lmmax*dx*np.sqrt(m2(lamb)), color='darkgray', alpha=0.05)
frmin = [np.argmin(np.abs(decay_times - lmmax))]
frmax = [np.argmin(np.abs(decay_times - lmmin))]
frmin = survive_prob(decay_times, ndcys)[frmin]
frmax = survive_prob(decay_times, ndcys)[frmax]
tmin = decay_times[np.argmin(np.abs(survive_prob(decay_times, ndcys) - frmax))]
tmax = decay_times[np.argmin(np.abs(survive_prob(decay_times, ndcys) - frmin))]
jfit_times = lin_fit_times(decay_times, ndcys, tmin, tmax)
gammas1[tmp] = np.array([temp, -jfit_times[0]])
col = next(cls)
ax[0].plot(decay_times*dx*np.sqrt(m2(lamb)), np.exp(get_line(decay_times, *jfit_times)), ls=':', lw=1, color='k')
ax[0].plot(decay_times*dx*np.sqrt(m2(lamb)), survive_prob(decay_times, ndcys), color=col, ls='-', label=labss)
temprange = np.linspace(tempList[0], tempList[3], 50)
bounce_dat0 = get_best_bounce(tempList[:4]/np.sqrt(m2(lamb)), gammas1[:4,1]/np.sqrt(m2(lamb)))
ax[1].plot(temprange/np.sqrt(m2(lamb)), exp_bounce(temprange/np.sqrt(m2(lamb)), *bounce_dat0), \
color='k', label=r'$A \sqrt{\frac{E}{T}} e^{-E/T}$')
temprange = np.linspace(save_teff_from_fld[0], save_teff_from_fld[3], 50)
bounce_dat = get_best_bounce(save_teff_from_fld[:4]/np.sqrt(m2(lamb)), gammas1[:4,1]/np.sqrt(m2(lamb)))
ax[2].plot(temprange/np.sqrt(m2(lamb)), exp_bounce(temprange/np.sqrt(m2(lamb)), *bounce_dat), \
color='k', label=r'$A \sqrt{\frac{E}{T_{\rm eff}}} e^{-E/T_{\rm eff}}$')
print(bounce_dat)
cls = cycle(allcolors[:4])
for tmp in range(len(tempList)):
if tmp>3: break
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
col = next(cls)
ax[1].plot(temp/np.sqrt(m2(lamb)), gammas1[tmp,1]/np.sqrt(m2(lamb)), color=col, marker='o', linestyle='none')
temp2 = save_teff_from_fld[tmp]
ax[2].plot(temp2/np.sqrt(m2(lamb)), gammas1[tmp,1]/np.sqrt(m2(lamb)), color=col, marker='o', linestyle='none')
ax[0].set_yscale('log')
#ax[0].ticklabel_format(axis='both', scilimits=[-5,5.])
ax[0].set_xlim(0, lmmax*dx*np.sqrt(m2(lamb)))
ax[0].set_ylim(4e-2, 1.05)
ax[0].set_ylabel(r'$\rm Pr(survive)$', fontsize=13)
ax[0].set_xlabel(r'$m \, t$')
h, l = ax[0].get_legend_handles_labels() # Extracting handles and labels
handles = [plt.plot([],marker="", ls="")[0]] + h
labels = [r'$T/m:$'] + l # Merging labels
leg = ax[0].legend(handles, labels, ncol=5, loc='lower left', labelspacing=0.3, columnspacing=0.6)#, markerfirst=False)
ax[0].grid(True, ls=':', color='lightgray', alpha=0.5)
ax[0].tick_params(direction='in', which='both', top=True, right=True)
for vpack in leg._legend_handle_box.get_children()[:1]:
for hpack in vpack.get_children():
for ii in hpack.get_children():
ii.set_width(-10)
for vpack in leg._legend_handle_box.get_children()[1:]:
for hpack in vpack.get_children():
for ii in hpack.get_children():
ii.set_width(18)
a = tempList[tmpList]/np.sqrt(m2(lamb))
ax[1].set_xticks(a)
a = [round(ii,2) for ii in save_teff_from_fld/np.sqrt(m2(lamb))]
ax[2].set_xticks(a)
for ai, aa in enumerate(ax[1:]):
aa.set_ylabel(r'$\Gamma / m \, L$', fontsize=13)
aa.set_xlabel([r'$T/m$',r'$T_{\rm eff}/m$'][ai])
aa.legend(loc=2)
aa.grid(True, ls=':', color='lightgray', alpha=0.5)
aa.tick_params(direction='in', which='both', top=True, right=True)
aa.ticklabel_format(axis='both', style='scientific', scilimits=[-1.,0.])
fig.tight_layout()
plt.savefig('./plots/decay_rate_and_surv_fraction.pdf')
plt.show()
[1.094875 3.09474548]
fig, ax = plt.subplots(1, 3, figsize = (9.7,3.), gridspec_kw={'width_ratios': [4, 2.3, 2.3]})
gammas1 = np.zeros((len(tempList), 2))
cls = cycle(allcolors[:4])
for tmp in range(len(tempList[:4])):
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
labss = r'${}$'.format(temp)
lmmax, lmmin = 5*nLat, 1.5*nLat#2*nLat/3
decay_times = np.load(decay_times_file(*exp_params, minSim, maxSim, nTimeMAX))
ndcys = maxSim - minSim
decay_times = np.sort(decay_times[:,1])
decay_times = decay_times[decay_times < lmmax]
ax[0].fill_betweenx(np.linspace(0,1.05,10), lmmin, lmmax, color='darkgray', alpha=0.05)
frmin = [np.argmin(np.abs(decay_times - lmmax))]
frmax = [np.argmin(np.abs(decay_times - lmmin))]
frmin = survive_prob(decay_times, ndcys)[frmin]
frmax = survive_prob(decay_times, ndcys)[frmax]
tmin = decay_times[np.argmin(np.abs(survive_prob(decay_times, ndcys) - frmax))]
tmax = decay_times[np.argmin(np.abs(survive_prob(decay_times, ndcys) - frmin))]
jfit_times = lin_fit_times(decay_times, ndcys, tmin, tmax)
gammas1[tmp] = np.array([temp, -jfit_times[0]])
col = next(cls)
ax[0].plot(decay_times, np.exp(get_line(decay_times, *jfit_times)), ls=':', lw=1, color='k')
ax[0].plot(decay_times, survive_prob(decay_times, ndcys), color=col, ls='-', label=labss)
temprange = np.linspace(tempList[0], tempList[3], 50)
bounce_dat0 = get_best_bounce(tempList[:4], gammas1[:4,1])
ax[1].plot(temprange, exp_bounce(temprange, *bounce_dat0), color='k', label=r'$\rm Best \; Fit$')#label=r'$A \sqrt{\frac{E}{T}} e^{-E/T}$')
temprange = np.linspace(save_teff_from_fld[0], save_teff_from_fld[3], 50)
bounce_dat = get_best_bounce(save_teff_from_fld[:4], gammas1[:4,1])
ax[2].plot(temprange, exp_bounce(temprange, *bounce_dat), color='k', label=r'$\rm Best \; Fit$')#label=r'$A \sqrt{\frac{E}{T_{\rm eff}}} e^{-E/T_{\rm eff}}$')
print(bounce_dat)
cls = cycle(allcolors[:4])
for tmp in range(len(tempList)):
if tmp>3: break
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
col = next(cls)
ax[1].plot(temp, gammas1[tmp,1], color=col, marker='o', linestyle='none')
temp2 = save_teff_from_fld[tmp]
ax[2].plot(temp2, gammas1[tmp,1], color=col, marker='o', linestyle='none')
ax[0].set_yscale('log')
ax[0].set_xlim(0, lmmax)
ax[0].set_ylabel(r'$\rm Pr(survive)$', fontsize=13)
ax[0].set_xlabel(r'$t$')
ax[0].grid(True, ls=':', color='lightgray', alpha=0.5)
ax[0].tick_params(direction='in', which='both', top=True, right=True)
a = tempList[tmpList]
ax[1].set_xticks(a)
a = [round(ii,3) for ii in save_teff_from_fld]
ax[2].set_xticks(a)
for ai, aa in enumerate(ax[1:]):
aa.set_ylabel(r'$\Gamma / L$', fontsize=13)
aa.set_xlabel([r'$T$',r'$T_{\rm eff}$'][ai])
aa.legend(loc=2)
aa.grid(True, ls=':', color='lightgray', alpha=0.5)
aa.tick_params(direction='in', which='both', top=True, right=True)
aa.ticklabel_format(axis='both', style='scientific', scilimits=[-1.,0.])
fig.tight_layout()
plt.show()
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in sqrt """Entry point for launching an IPython kernel. /cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: overflow encountered in exp """Entry point for launching an IPython kernel.
[0.10948856 0.30947499]
def get_stuffs(nsets):
partdat = np.zeros((len(tmpList), nsets))
for ii, average_bubble in enumerate(bubbleList):
tmp = tmpList[ii]
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
veldata = np.array(np.load(velocities_file(*exp_params)))
simvels, all_vels = veldata[:,0].astype(int), veldata[:,1]
pool = set(all_vels)
means, vars = np.empty(nsets), np.empty(nsets)
slen = int(round(len(pool) / nsets,1))
for ss in range(nsets):
velssec = random.sample(pool, slen)
pool -= set(velssec)
velssec = np.array(velssec)
partdat[tmp,ss] = np.std(velssec)
return partdat
ecrit, esph, char_vels = np.zeros((3, len(tmpList)))
for ii, average_bubble in enumerate(reversed(bubbleList)):
tmp = tmpList[ii]
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
veldata = np.array(np.load(velocities_file(*exp_params)))
simvels, all_vels = veldata[:,0], veldata[:,1]
char_vels[tmp] = np.std(all_vels)
instanton = instantonList[tmp]
gradinst = (np.roll(instanton,-1) - instanton)/dx
ecrit[tmp] = np.sum(0.5*gradinst**2. + V(instanton, lamb) - V(phieq,lamb))
esph[tmp] = encritList[tmp]
diffen = phi0**2.*dx
#sigma_avbub = tempList[tmpList] / (esph * diffen)
#sigma_crit = tempList[tmpList] / (ecrit * diffen)
#sigma_drate = tempList[tmpList] / bounce_dat[-1]
sigma_avbub = save_teff_from_fld / (esph * diffen)
sigma_crit = save_teff_from_fld / (ecrit * diffen)
sigma_drate = save_teff_from_fld / bounce_dat[-1]
fig, ax = plt.subplots(1,4, figsize = (8,2), sharey='row')
cols = cycle(allcolors[:4])
save_gauss_fit = []
for ii, average_bubble in enumerate(bubbleList):
tmp = tmpList[ii]
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
labss = r'$T/m = {}$'.format(round(temp/np.sqrt(m2(lamb)),1))#labl(lamb, phi0, temp)
veldata = np.array(np.load(velocities_file(*exp_params)))
simvels, all_vels = veldata[:,0].astype(int), veldata[:,1]
col = next(cols)
n, bins, patches = ax[ii].hist(all_vels, bins=20, alpha=1, color=col, density=True)#, label=labss)
n = n[1:]
bins = bins[1:]
centers = 0.5*(bins[1:] + bins[:-1])
xlist = np.linspace(-1.,1.,1000)
gauss_fit = lambda x, mean, std: np.exp(- (x - mean)**2. / (2.*std**2.)) / np.sqrt(2.*np.pi*std**2.)
fit_gauss = lambda x, data: sco.curve_fit(gauss_fit, x, data)[0]
best_gauss = fit_gauss(centers, n)
ax[ii].plot(xlist, gauss_fit(xlist, *best_gauss), color='k', ls='--')
save_gauss_fit.append(best_gauss[-1])
bb = sigma_crit[tmp]**0.5
ax[ii].plot(xlist, gauss_fit(xlist, 0., bb), color='k', ls='-')
ax[ii].plot(xlist, gauss_fit(xlist, np.mean(all_vels), np.std(all_vels)), color='k', ls=':')
a = ax[ii].get_xticks().tolist()[1:-1:]
a = np.linspace(-1, 1, 5)
a[::2] = np.array(a[::2], dtype='int')
ax[ii].set_xticks(a)
b = [r'${:.1f}$'.format(al) for al in a]
b[::2] = [r'${:.0f}$'.format(al) for al in a[::2]]
ax[ii].set_xticklabels(b)
ax[ii].grid(ls=':', color='darkgray',alpha=0.5)
ax[ii].set_xlabel(r'$v$')
ax[ii].set_title(labss, fontsize=10)
ax[ii].tick_params(direction='in', which='both', top=True, right=True)
a = ax[0].get_yticks().tolist()[:-1:]
a = [round(al,2) for al in a]
ax[0].set_yticks(a)
b = [r'${:.1f}$'.format(al) for al in a]
ax[0].set_yticklabels(b)
ax[0].set_ylabel(r'$\rm PDF$')
plt.savefig('./plots/comparison_vels_PDF.pdf')
plt.show()
fig, ax = plt.subplots(1,1, figsize = (4,2.5))
ax.plot(save_teff_from_fld/np.sqrt(m2(lamb)), sigma_crit, label=r'$\rm Prediction$')
ax.plot(save_teff_from_fld/np.sqrt(m2(lamb)), sigma_avbub, label=r'$E_{\rm exp}$')
ax.plot(save_teff_from_fld/np.sqrt(m2(lamb)), sigma_drate, label=r'$\rm Decay \; Rate$')
for nsets in reversed(np.linspace(15, 20, 1)):
partdat = get_stuffs(int(nsets))
ax.errorbar(save_teff_from_fld/np.sqrt(m2(lamb)), np.mean(partdat, axis=-1)**2., yerr=np.std(partdat, axis=-1), \
alpha=1, lw=1, label=nsets)
ax.plot(save_teff_from_fld/np.sqrt(m2(lamb)), char_vels**2., label=r'$\rm Data$')
a = [round(ii,3) for ii in save_teff_from_fld/np.sqrt(m2(lamb))]
ax.set_xticks(a)
ax.grid(True, ls=':', color='lightgray', alpha=0.5)
ax.tick_params(direction='in', which='both', top=True, right=True)
ax.ticklabel_format(axis='both', style='scientific', scilimits=[-1.,0.])
ax.legend(bbox_to_anchor=(1, 1))
ax.set_xlabel(r'$T_{\rm eff}/m$')
ax.set_ylabel(r'$\left<v^2\right>$')
plt.tight_layout()
plt.show()
fig, ax = plt.subplots(1,1, figsize = (5.5,3))
ax.plot(save_teff_from_fld/np.sqrt(m2(lamb)), sigma_crit, label=r'$\rm Prediction$')
partdat = get_stuffs(15)
ax.errorbar(save_teff_from_fld/np.sqrt(m2(lamb)), char_vels**2., \
yerr=np.std(partdat, axis=-1), alpha=1, lw=1, label=r'$\rm Variance \; of \; Data$')
ax.errorbar(save_teff_from_fld/np.sqrt(m2(lamb)), np.array(save_gauss_fit)**2., \
yerr=np.std(partdat, axis=-1), alpha=1, lw=1, label=r'$\rm Gaussian \; Best \; Fit$')
a = [round(ii,3) for ii in save_teff_from_fld/np.sqrt(m2(lamb))]
ax.set_xticks(a)
ax.grid(True, ls=':', color='lightgray', alpha=0.5)
ax.tick_params(direction='in', which='both', top=True, right=True)
ax.ticklabel_format(axis='x', style='scientific', scilimits=[-1.,0.])
ax.ticklabel_format(axis='y', style='scientific', scilimits=[-10.,10.])
a = ax.get_yticks().tolist()[1::2]
a = [round(al,2) for al in a]
ax.set_yticks(a)
a = [r'${:.2f}$'.format(al) for al in a]
ax.set_yticklabels(a)
leg = ax.legend(bbox_to_anchor=(1.5, 0.6))._legend_box.align='left'
ax.set_xlabel(r'$T_{\rm eff}/m$')
ax.set_ylabel(r'$\left<v^2\right>$')
plt.tight_layout()
plt.savefig('./plots/comparison_vcom_pred.pdf')
plt.show()
#velsq_exp = lambda edt: 1. - edt - edt**2.* np.exp(edt) * scp.special.expi(-edt)
#velsq_exp_explicit = lambda en, tm: 1. - en/tm - (en/tm)**2.* np.exp(en/tm) * scp.special.expi(-en/tm)
#def f_inv_velsq_exp_explicit(x, ii):
# return abs(velsq_exp_explicit(x, tempList[tmpList][ii]) - char_vels[ii])
#enfromvels = np.zeros((len(char_vels)))
#for ii in range(len(char_vels)):
# out = sco.minimize_scalar(f_inv_velsq_exp_explicit, bounds=((0,None)), args=(ii))
# enfromvels[ii] = out.x
# Plot to compare all energies
fig, ax = plt.subplots(1, 1, figsize = (5.5,3.))
ax.plot(save_teff_from_fld/np.sqrt(m2(lamb)), 1./sigma_crit, label=r'$\rm Prediction$')
ax.plot(save_teff_from_fld/np.sqrt(m2(lamb)), 1./sigma_avbub, label=r'$\rm Critical \; Bubble$')
ax.plot(save_teff_from_fld/np.sqrt(m2(lamb)), 1./sigma_drate, label=r'$\rm Decay \; Rate$')
ax.plot(save_teff_from_fld/np.sqrt(m2(lamb)), 1./np.array(save_gauss_fit)**2., label=r'$\rm Velocity \; Distribution$')
a = [round(ii,3) for ii in save_teff_from_fld/np.sqrt(m2(lamb))]
ax.set_xticks(a)
ax.ticklabel_format(axis='x', style='scientific', scilimits=[-1.,0.])
ax.ticklabel_format(axis='y', style='scientific', scilimits=[-10.,10.])
ax.set_ylabel(r'$E/T_{\rm eff}$')
ax.set_xlabel(r'$T_{\rm eff}/m$')
ax.legend(ncol=1, loc='best', frameon=False, bbox_to_anchor=(1, 1))._legend_box.align='left'
ax.grid(which='both', color='darkgray', ls=':', alpha=0.4)
ax.tick_params(direction='in', which='both', top=True, right=True)
plt.tight_layout()
plt.savefig('./plots/critenergy_comparison.pdf')
plt.show()
for ii, average_bubble in enumerate(bubbleList):
tmp = tmpList[ii]
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
decay_times = np.load(decay_times_file(*exp_params, minSim, maxSim, nTimeMAX))
print('N1 = ', len(decay_times))
minDecTime = nLat*2//3
alltimes = decay_times[:,1]
simList2Do = decay_times[alltimes>=minDecTime, 0]
n2Do = len(simList2Do)
print('N2 = ', n2Do)
veldata = np.array(np.load(velocities_file(*exp_params)))
simvels, all_vels = veldata[:,0].astype(int), veldata[:,1]
print('N3 = ', len(all_vels))
# print(temp, len(all_vels), np.count_nonzero(all_vels > 0.)/len(all_vels))
# print(temp, np.count_nonzero(np.abs(all_vels) > 0.5)/len(all_vels))
N1 = 1003 N2 = 665 N3 = 654 N1 = 2192 N2 = 1199 N3 = 1191 N1 = 3165 N2 = 1280 N3 = 1266 N1 = 3711 N2 = 882 N3 = 854
if False:
for tmp in range(len(tempList)):
if tmp!=0: continue
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
# 0 for field; 1 for momentum
find = 0
# modes to plot:
aa, bb = 1, knyq-1
labs = labl(lamb, phi0, temp)
psth = pspec(lamb, phi0, temp)[aa:bb]
ALL_powspec1 = np.load(powspec_tlist_file(*exp_params, minSim, 2000))
ALL_powspec2 = np.load(powspec_tlist_file(*exp_params, 2000, maxSim))
tlist, PSfld1, PSfld2 = ALL_powspec1[0], ALL_powspec1[1], ALL_powspec2[1]
fig, ax = plt.subplots(1,1, figsize = (8,4))
camera = Camera(fig)
for tind, tt in enumerate(tlist):
avPSfld1 = np.nanmean(PSfld1[:, find, tind, :], axis=0)
avPSfld2 = np.nanmean(PSfld2[:, find, tind, :], axis=0)
avPSfld = np.nanmean(np.array([avPSfld1[:], avPSfld2[:]]), axis=0)
slice = plt.plot(klist[aa:bb], avPSfld[aa:bb]/psth, ls='-', color='purple')
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_ylabel(r'$P(k)$')
ax.set_xlabel(r'$k$')
ax.legend(slice, [f't = {tt}'], loc=1, title=labs)
camera.snap()
animation = camera.animate(interval = 0.05);
animation.save('./plots/animation_cut_PS'+batch_params(*exp_params)+'.gif', writer = 'imagemagick')
fig, ax = plt.subplots(1,1, figsize = (5,2.5))
for tmp in range(len(tempList)):
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
aa, bb = 1, knyq-1
plt.plot(klist[aa:bb], pspec(lamb, phi0, temp)[aa:bb], ls=lsl(tmp), label = labl(lamb, phi0, temp)) # th pow spec
plt.xscale('log')
plt.legend(); plt.show()
if True:
choose = random.sample(np.arange(4000).tolist(), 4)
for sim in choose:
fig, ax = plt.subplots(1,1, figsize = (4,2.5))
for tmp in range(len(tempList)):
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
if sim < 2000: minSim, maxSim = 0, 2000
else: minSim, maxSim, sim = 2000, 4000, sim-2000
ALL_toten = np.load(toten_tlist_file(*exp_params, minSim, maxSim))
tlist, energy = ALL_toten[0], ALL_toten[1][sim]
try: nnrg = np.argwhere(np.isnan(energy)).flatten()[0]
except: nnrg = len(tlist)+1
tcut, tencut = tlist[:nnrg], energy[:nnrg]
tencut = dx * phi0**2. * (tencut - tencut[0]) / tencut[0]
plt.plot(tcut, tencut, label=labl(lamb, phi0, temp))
plt.legend(title='Sim='+str(sim), bbox_to_anchor=(1.2,1.))
plt.show()
if True:
for tmp in range(len(tempList)):
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
labs = labl(lamb, phi0, temp)
minSim = 2000
ALL_toten = np.load(toten_tlist_file(*exp_params, minSim, maxSim))
ALL_emt = np.load(emt_tlist_file(*exp_params, minSim, maxSim))
tlist, enfld, emtfld = ALL_toten[0], ALL_toten[1], ALL_emt[1]
fig, ax = plt.subplots(1,1, figsize = (4,2.5))
for sim, (momentum, energy) in enumerate(zip(emtfld, enfld)):
if sim!=0: continue
try: nnrg = np.argwhere(np.isnan(energy)).flatten()[0]
except: nnrg = len(tlist)+1
tcut, emtcut, tencut = tlist[:nnrg], momentum[:nnrg], energy[:nnrg]
emtcut = (emtcut - emtcut[0]) / emtcut[0]
tencut = (tencut - tencut[0]) / tencut[0]
ax.plot(tcut, tencut, label=(r'$T^{00}$' if sim==0 else None))
ax.plot(tcut, emtcut, label=(r'$T^{0x}$' if sim==0 else None))
ax.set_xlabel(r'$\phi_0^{-1} \sqrt{V_0} \; t$')
ax.set_ylabel(r'$T^{00}, \; T^{0x}$')
ax.legend(title=labs, bbox_to_anchor=(1.,1.))
ax.grid(ls=':', color='darkgray', alpha=0.3)
ax.set_title(labs)
fig.show()
if True:
fig, ax = plt.subplots(1, 1, figsize=(3.5, 2.5))
allcs = cycle(allcolors)
for tmp in reversed(range(len(tempList))):
if tmp>3: continue
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
labss = r'${}$'.format(round(temp/np.sqrt(m2(lamb)),1))#labl(lamb, phi0, temp)
col = next(allcs)
veldata = np.array(np.load(velocities_file(*exp_params)))
simvels, all_vels = veldata[:,0].astype(int), veldata[:,1]
ALL_emt1 = np.load(emt_tlist_file(*exp_params, minSim, 2000))
ALL_emt2 = np.load(emt_tlist_file(*exp_params, 2000, maxSim))
emtfld = np.concatenate((ALL_emt1[1], ALL_emt2[1]), axis=0)
initemt = emtfld[simvels,0]
#emtfld = np.load(stdemt0_tlist_file(*exp_params, minSim, maxSim))
#initemt = emtfld[simvels]
ax.plot(all_vels, -initemt, marker='o', ms=2, linestyle='None', label=labss, color=col)
f = lambda x, A: x/A
pbf, pcov = sco.curve_fit(f, all_vels, -initemt) # your data x, y to fit
ax.plot(all_vels, f(all_vels, *pbf), color='k', linewidth=1)
ax.set_xlim(-1,1)
# ax.set_ylim(-1,1)
ax.set_xlabel(r'$v_{\rm COM}$')
# ax.set_ylabel(r'$\int \mathrm{d} r \, \Pi \, \partial_r \varphi$')
ax.set_ylabel(r'$\bar{P}$')
ax.grid(ls=':', color='darkgray', alpha=0.3)
leg = ax.legend(title=r'$T/m:$', ncol=1, loc='best', fontsize=10, bbox_to_anchor=(1, 0.8))
if True:
a = ax.get_yticks().tolist()[2:-1:2]
a = [round(al,2) for al in a]
ax.set_yticks(a)
# a[::2] = [int(al) for al in a[::2]]
# b = [r'${:.1f}$'.format(al) for al in a]
# b[::2] = [r'${:.0f}$'.format(al) for al in a[::2]]
# ax.set_yticklabels(b)
a = ax.get_xticks().tolist()
a = [round(al,2) for al in a]
ax.set_xticks(a)
a[::2] = [int(al) for al in a[::2]]
b = [r'${:.1f}$'.format(al) for al in a]
b[::2] = [r'${:.0f}$'.format(al) for al in a[::2]]
ax.set_xticklabels(b)
plt.savefig('./plots/initial_emt_vs_vel.pdf')
fig.show()
if True:
fig, ax = plt.subplots(1, 1, figsize=(5, 4))
allcs = cycle(allcolors)
for tmp in reversed(range(len(tempList))):
if tmp>3: continue
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
labss = r'${}$'.format(round(temp/np.sqrt(m2(lamb)),1))#labl(lamb, phi0, temp)
col = next(allcs)
veldata = np.array(np.load(velocities_file(*exp_params)))
simvels, all_vels = veldata[:,0].astype(int), veldata[:,1]
ALL_emt1 = np.load(emt_tlist_file(*exp_params, minSim, 2000))
ALL_emt2 = np.load(emt_tlist_file(*exp_params, 2000, maxSim))
emtfld = np.concatenate((ALL_emt1[1], ALL_emt2[1]), axis=0)
initemt = emtfld[simvels,0]
# fit a straight line
jfit_times = np.polyfit(all_vels, -initemt, deg=1)
print(jfit_times)
vlist = np.linspace(0., 1, 100)
ax.plot(vlist, get_line(vlist, *jfit_times), '-', label=labss, color=col)
ax.legend(title=r'$T/m {\rm \; best \; fit}$', loc='best', ncol=2)
ax.set_xlabel(r'$v_{\rm COM}$')
ax.set_ylabel(r'$\bar{P}$')
ax.grid(ls=':', color='darkgray', alpha=0.3)
fig.show()
[ 0.18583269 -0.00668639] [ 0.18248509 -0.00590275] [0.18450292 0.00386879] [0.18669835 0.00104476]
if False:
# for tmp in reversed(range(len(tempList))):
for tmp in reversed(range(len(tempList))):
if tmp>3: continue
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
path = decay_times_file(*exp_params, minSim, maxSim, nTimeMAX)
if os.path.exists(path):
print(path)
decay_times = np.load(path)
minDecTime = nLat*2//3
alltimes = decay_times[:,1]
simList2Do = decay_times[alltimes>=minDecTime, 0]
all_vels = []
for sim in simList2Do:
path2RESTsim = rest_sim_location(*exp_params, sim)
if os.path.exists(path2RESTsim):
sim, bubble, totbeta = np.load(path2RESTsim)
all_vels.append(np.array([sim, totbeta]))
np.save(velocities_file(*exp_params), all_vels)
print(len(all_vels), len(simList2Do))
print('Done!')
if True:
fig2, ax = plt.subplots(1, 1, figsize=(6, 5))
plt.grid(True, ls=':', alpha=0.5)
for tmp in reversed(range(len(tempList))):
if tmp>3: continue
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
exp_params = np.asarray([nLat, lamb, phi0, temp])
labss = r'${}$'.format(round(temp/np.sqrt(m2(lamb)),1))#labl(lamb, phi0, temp)
veldata = np.load(velocities_file(*exp_params))
simvels, all_vels = veldata[:,0].astype(int), veldata[:,1]
print(temp, len(simvels))
plt.hist(all_vels, bins=12, label=labss, density=True)
plt.xlim((-1,1))
plt.legend(title=r'$T/m$', ncol=1, loc='best')
plt.xlabel(r'$v$')
plt.legend()
plt.show()
0.12 854 0.11 1266 0.1 1191 0.09 654